arXiv:1505.03250vl [math.NA] 13 May 2015 


Multiscale numerical schemes for kinetic equations in the 

anomalous diffusion limit 


Nicolas Crouseilles a,b , Helene Hivert b Mohammed Lemou c,b 

a INRIA 

b IRMAR, Universite de Rennes 1. Campus de Beaulieu. 35000 Rennes. 

C CNRS 

Received *****; accepted after revision +++++ 

Presented by 


Abstract 


We construct numerical schemes to solve kinetic equations with anomalous diffusion scaling. When the equilib¬ 
rium is heavy-tailed or when the collision frequency degenerates for small velocities, an appropriate scaling should 
be made and the limit model is the so-called anomalous or fractional diffusion model. Our first scheme is based on 
a suitable micro-macro decomposition of the distribution function whereas our second scheme relies on a Duhamel 
formulation of the kinetic equation. Both are Asymptotic Preserving (AP): they are consistent with the kinetic 
equation for all fixed value of the scaling parameter e > 0 and degenerate into a consistent scheme solving the 
asymptotic model when e tends to 0. The second scheme enjoys the stronger property of being uniformly accurate 
(UA) with respect to e. The usual AP schemes known for the classical diffusion limit cannot be directly applied 
to the context of anomalous diffusion scaling, since they are not able to capture the important effects of large 
and small velocities. We present numerical tests to highlight the efficiency of our schemes. To cite this article: N. 
Crouseilles, H. Hivert, M. Lemou, C. R. Acad. Sci. Paris, Ser. I 3^0 (2005) ???????. 

Resume 

Schemas numeriques multi-echelles pour les equations cinetiques dans la limite de diffusion anor- 
male. 

Nous construisons des schemas numeriques pour resoudre les equations cinetiques dans le regime de diffusion 
anormale. Lorsque l’equilibre presente une queue lourde ou lorsque la frequence de collision degenere pour les 
petites vitesses, un scaling approprie permet d’obtenir un modele asymptotique appele modele de diffusion anor¬ 
male ou fractionnaire. Le premier schema que nous construisons est base sur une decomposition micro-macro de 
la fonction de distribution tandis que le second s’appuie sur une formulation de Duhamel de l’equation de depart. 
Ces deux schemas sont Asymptotic Preserving (AP) : ils sont consistants avec l’equation cinetique lorsque le pa- 
rametre d’echelle e > 0 est fixe et degenerent en un schema consistant avec le modele limite quand e tend vers 0. 
Le deuxieme schema est meme uniformement precis (UA) par rapport a e. Les schemas AP qui sont connus dans 
le cas de la limite de diffusion classique ne peuvent pas directement s’appliquer au cas de la diffusion anormale 
car ils ne permettent de capturer les effets importants des petites et des grandes vitesses. Nous presentons des 
tests numeriques pour mettre en evidence l’efficacite des schemas que nous presentons. 
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Version franqaise abregee 

Le but de ce travail est de mettre en place des schemas numeriques pour les equations cinetiques 
lineaires dans le cas de la limite de diffusion anonrale. Quand la distribution d’equilibre est une fonction 
a queue lourde ou lorsque la frequence de collision est degeneree pour les petites vitesses, l’equation 
cinetique (1) tend vers l’equation dite de diffusion anormale (2) lorsque le parametre d’echelle e tend 
vers zero. Dans cette limite, une raideur apparait dans l’equation de depart et une resolution numerique 
directe du probleme peut devenir extremement couteuse puisque les parametres numeriques doivent a 
priori etre adaptes a e. La construction de schemas dits Asymptotic Preserving (AP) permet de repondre 
a cette contrainte : ces schemas restent consistants avec l’equation cinetique tout en s’affranchissant de 
la contrainte sur les parametres numeriques, et degenerent vers l’equation limite quand e tend vers 0. 

Dans le cas de l’asymptotique de diffusion anormale, la rnise en place de tels schemas s’avere plus 
compliquee que dans le cas classique. En effet, en plus de la raideur evoquee ci-dessus (e tend vers 0), il 
est crucial de capturer les effets des grandes et petites vitesses pour que ces schemas degenerent vers des 
approximations consistantes de l’equation de diffusion anormale (2). Un schema numerique inspire des 
approclies AP standards ne tiendrait pas compte de ces effets, et degenererait vers une approximation 
d’une equation de diffusion classique et non vers celle du rnodele correct de diffusion anormale. 

Deux cas sont consideres dans ce travail : le cas d’un equilibre a queue lourde et celui d’une frequence 
de collision degeneree en 0. Dans les deux cas, nous construisons d’abord un schema base sur une 
decomposition micro-macro de la solution, qui fournit un schema multi-echelle pour (1) completement ex- 
plicite en temps; ensuite, nous presentons un schema base sur une formulation de Duhanrel de l’equation 
cinetique (1) ayant une propriety plus forte que la propriete AP : la precision de ce schema est uniforme 
(UA) par rapport a e. 

Dans chaque cas, l’ecriture d’un schema semi-discret en temps permet d’obtenir une formulation qui 
tend vers l’equation de diffusion anormale lorsque e tend vers 0 si l’espace des vitesses est considere conune 
continu et les integrations en vitesses sont realisees exactenrent. Cependant, une discretisation directe de 
l’espace des vitesses pour realiser les integrations numeriquement ne permet pas de prendre en compte 
les effets des grandes et petites vitesses, qui sont a l’origine de la limite de diffusion anormale. Dans ce 
travail, nous montrons qu’il est done necessaire d’effectuer des transformations sur certaines integrates en 
vitesses avant de les discretiser. En l’occurrence, nous effectuons des changements de variables adequats 
pour faire apparaitre naturellement les termes a l’origine de l’equation asymptotique. Nous obtenons ainsi 
des schemas completement discretises ayant la propriete AP et UA. Nous presentons egalement des tests 
numeriques qui mettent en evidence la limite de diffusion anormale de nos schemas quand £ tend vers 0. 
Cette note est une version abregee de [3], [4]. 


1. Introduction 

We consider the kinetic equation 


Email addresses: nicolas.crouseilles@inria.fr (Nicolas Crouseilles), helene.hivert@univ-rennesl.fr (Helene 
Hivert), mohanraied.lemou@univ-rennesl.fr (Mohammed Lemou). 
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£ a dtf + £V V x f = Q(f) := i/(v) {p v M -/), f(0,x,v) = f o (x,v), (1) 

where / is a distribution function which depends on the time t > 0, the space variable x £ and the 
velocity v £ with d = 1,2,3, fo is a given initial data. In the sequel, we will denote by brackets the 
integration in v. We define the density p by p(t, x) = ( f(t , x, v)) =: f v£R d f(t, x, v)dv, and the quantity p u 
(which is not the usual density) by p v {t, x) = ( v{v)f) / {v(v)M{v)). Note that this definition of p v ensures 
the local mass conservation ( Q(f )) = 0. The positive number e is a scaling parameter and a is a power 
which will be chosen according to the physical nature of the problem (see below for more details). The 
equilibrium distribution function M is a normalized even function of v. We will consider two physically 
relevant cases, depending on the nature of the equilibrium function M (see [6], [2]) and of the collision 
frequency v (see [1]): 

(i) Case 1 - heavy-tail: M is a power-tailed equilibrium and v(y) = 1. To simplify, we will consider 
the case M(v) = m/{ 1 + \v\@) where /3 £ (d, d + 2) and the normalization parameter m is chosen 
such that (M) = 1. In this case, the appropriate choice of a is: a = /3 — d G (0,2). 

(ii) Case 2 - degenerate collision frequency: v(y) ~ vo\v\ d+2+l3 , (3 > 0, > 0 and M(v) = 

v —^0 

e ^2 / ' simplify^ we will consider the case v(v) = t'o|^| d+2+/3 for some (5 > 0. In this case, the 
appropriate choice of a is: a= (2 + 2d + /3)/(l + d + f3) £ (1, 2). 

Note that in these two cases, the quantity D = which usually appears as the diffusion coefficient 

in the classical diffusion case, is not finite. In fact, in our context, the limit for small £ of (1) is not given 
by diffusion but by anomalous diffusion equation, which can be written with a fractional Laplacian 

d t p + K \k\ a p = 0, p(0,k) = ^/o(M)) • (2) 

In (2), the quantity p (respectively /) denotes the space Fourier transform of the function p (respectively 
/) and k is the Fourier variable. The coefficient k can be expressed in the two cases: in Case 1 of heavy-tail 
equilibrium, it writes 

/to (v ■ e) 2 \ 

\ \v\P 1 + (v ■ e) 2 / ’ 

and in Case 2 of the degenerate collision frequency, it writes 

= 1 ~ a ( 1 {v ■ e) 2 \ 

v^d+l + /3 \ \v\P l + ive) 2 / 1 

where e denotes any unitary vector of R d . The anomalous diffusion limit for kinetic equations has been 
studied in [1] in the case of degenerate collision frequency and in [2], [6] in the case of heavy-tailed 
equilibrium. 


2. Micro-macro scheme 

In this section, we derive an AP micro-macro scheme for (1) in the case of the anomalous diffusion limit. 
This approach bears similarities with the one developed in [5] and has the strong advantage of treating 
the transport explicitly. Note that an implicit AP scheme can be constructed directly, with no use of the 
micro-macro decomposition. Indeed we proceed as follows. We first start by an implicit formulation of 
(1) 

r +1 = A(u)(^/+^^V^ pZ +1 M + (i - \{v)) (V + ^yv ■ (3) 
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where X(v) = Atv(v)/(£ a + A tv(v)), At is the time step, f n ~ f{t n ), p n = (/”) and g n = f n — p n M such 
that (c/ 11 ) = 0. We then integrate (3) against v with respect to v and get an expression of p” +1 ; then, we 
use an adequate treatment of the integrals in velocity appearing in the so-obtained expression of p” +1 
to ensure the AP property of the scheme, using adapted changes of variables (see [3] for more details). 
Once this expression of p" +1 is obtained, we report it in (3) to compute / n+1 . Here, we are interested in 
deriving numerical schemes where the transport part is explicit in time. To this end, we use a suitable 
micro-macro decomposition which leads to the following numerical scheme (see again [3] for more details). 
Proposition 2.1 We introduce the semi-discrete micro-macro scheme defined for all x € K d ,i> £ and 
all time index 0 < n < N, with NAt = T (T > 0), by 



n 



+ £ 


f'V/ \ 
e a + v{v)At / 


= 0, 


,»+1 


At 


+ E 1 -*v ■ V x p n M + £ x - a (V ■ V x g n - (v ■ S7 x g n ) M) = 


^ (»” +1 


(4) 


{v{v)g n+1 ) 
(v{v)M) 



(5) 


where J-~ l denotes the inverse Fourier transform in space. The quantity p n+1 / 2 can he chosen equal to 
p n or to p n+1 depending on the desired asymptotic scheme (explicit or implicit in time) and the quantity 
I is given by 


I = £ 


1—cx 


\{v)- 


i k ■ v 


1 + 


£ X(v] 

~v(v) 


i k ■ v 


-M ) = £ 


u(v)X(v) 2 (k ■ v ) 2 
v{v) 2 + £ 2 X(v) 2 (k ■ v ) 2 


M 


( 6 ) 


This scheme is of order 1 for any fixed £ > 0 and enjoys the AP property: for a fixed At, the scheme 
degenerates into a first order in time scheme for (2) when £ goes to zero. 

Let us say a few words about the derivation of (4)-(5). Since A = O(At), we have from (3) 


r +i = A(u) /+ 


eX(v) 

v(v) 


vV x ) p^M + (l-X(v))f n + 0(At). 


(7) 


Then, we integrate (1) in v to get the continuity equation on p and write a Euler implicit scheme: 
p n+1 = p n — Ate 1_a V x • (vf n+1 ). Then, we replace f n+1 by (7) in this scheme, use the identity p" +1 = 
pn+i (yg n+1 )I(yM) and the evenness of M(v) to get (4). To get (5), we just apply (I— n) to (1) where 
n/ = ( f)M and discretize the obtained equation on g by an explicit scheme for the transport terms and 
an implicit scheme for the collision terms. 

Before discussing the delicate issue relative to the velocity discretization in (4)-(5), we first briefly 
explain how ( p n ,g n ) are computed recursively from (4)-(5), assuming that the space and velocity dis¬ 
cretizations have been already fixed. The idea is to start with (5) to find an expression for g n+1 . In Case 1 
(heavy-tailed equilibrium) as v(v) = 1, the term (y(v)g n+1 ) in the right hand side of (5) vanishes, giving 
easily an expression for g n+1 , which is then reported in (4) and so on. However, in Case 2 (degenerate 
collision frequency), it is necessary to extract (y(v)g n+1 ) before solving the equation in g n+1 . To do that, 
we express g n+l from (5) in terms of p n ,g n and {y(v)g n+l )\ we multiply this obtained expression of g n+1 
by v(v) and integrate in v to get the following expression for (y(v)g n+1 ^ 


Hv)g n+1 ) 


/ Hv)g n \ / ^v)(^vV x p n M + ^(v.X/ x g n ^(y.V x g n )M)) 

\l + J±v(v) / \ 1 + 

At / v(v) 2 M \ 

£ a ( v(y)M ) \ 1 + j£-v{v) / 


Reporting this in (5) enables to compute g n+1 , which is reported into (4) to get p n+1 . 
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Now, we come to the construction of complete discretization from (4)-(5) and in particular the important 
point of discrete integrations in velocity. Indeed, the behavior of the scheme when e goes to 0 is intimately 
related to the discretization we make to approximate (6). To see this, we observe that a direct discretization 
(by rectangle formula for instance) of I given by (6) converges to 0 when e goes to 0 since the numerical 
approximation of the bracket in (6) is always finite. Therefore this leads to the wrong limit and to overcome 
this problem, we proceed as follows. The general idea is to perform a suitable change of variables in (6) 
before discretizing it in velocity. In Case 1 (heavy-tailed equilibrium) we make the change of variable 
w = e\\k\v in I before applying the discretization and in Case 2 (degenerate collision frequency) we make 
the change of variables w = e\k\v/v{v) in I before discretizing the brackets. Hence, when e tends to 0 ,1 
degenerates into the discretized coefficient k of (2). This ensures the AP property of the fully discretized 
scheme. Note that the other velocity integrations in (4)-(5) are discretized directly without any change 
of variables. Finally, a standard upwind scheme is used for the spatial discretization of (5). 

To highlight the AP character of this scheme, we present in Fig. 1 the densities p computed at time 
T = 0.1 with this scheme for some e. We took the initial data fo{x, v ) = (1 + sin(7 tx))M(v) for x G [—1,1] 
and 64 points of discretization in space. In Case 1, we took /? = 2.5 and in Case 2 we took /? = 0.5. In 
both cases we considered u max = 5, and 200 points of discretization in velocity and At = 10 -3 . 




Figure 1. The densities p(T = 0.1. x) for the micro-macro schemes and the anomalous diffusion limit scheme. Left: the case 
of heavy-tailed equilibrium. Right: the case of degenerate collision frequency. 


3. Duhamel formulation based scheme 

In this section, our goal is to derive numerical schemes for (1) which go beyond the AP property. More 
precisely, our aim to construct numerical schemes whose accuracy in time is uniform with respect to e. 
Note that although the scheme of the last section is AP, it does not enjoy this uniform accuracy property 
(see [3]). The starting point of the derivation is a Duhamel formulation of (1) in Fourier variable, which 
leads, after an integration in v against is, to an exact equation on p v 

t 

p v {t, k) = A 0 (t, k ) + ((is(v)A'I))~ 1 j i s[v) 2 M(v)^ p v (t - s a s, k)ds, (8) 

with A 0 (t,k ) = (e~'P r< ' u ^ +1£k " v )is(v)fo(k, V )'^ / (is(v)M). Considering discrete time values t n = nAt,0 < 
n < N, NAt = T, ( T , At > 0), and evaluating (8) at t = t n+ \ leads to 


For s G [tj/e a , tj+i/e a ], using p” ~ p v {t n ), we consider the approximation p v {t n +\ — e a s,k ) « (tj+i — 
£ a s)/At + (e“s — tj)/At p™- j and get the following proposition. 


p v {t n +i,k) = A 0 (t n+ i,k) + (( is(v)M)) 1 y, 

i=o 


J tj e (e- s ^vW-v) v (vYM{v)) pAtn+l - £ a s, k)ds. 
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Proposition 3.1 We introduce the Duhamel formulation based scheme defined for all x £ £ R. d 

and all time index 0 < n < N, NAt = T by 


A n +1 


= A 0 (t n+1 ,k) + 


_-_V 

{u(v)M) 


( c j - b j)P, 


n+l-j 


b jP 


n~3 


(9) 


where 


*3 +1 *3 + 1 



This scheme is of order 1 uniformly in e: 3 C > 0 independent from e, such that maxj |/3„(f™, k ) — pf,(k)\ < 
CAt,\/0<n<N,NAt = T. 

We refer to [3] and [4] for more details and for the proof of the uniform accuracy of this scheme. 

The result of this proposition implies in particular, that this semi-discrete scheme enjoys the AP 
property. Moreover, we will show below how to discretize the velocity integrals in order to preserve this 
uniform accuracy property for the fully discretized scheme. In fact, as in the previous section, a direct 
rectangular velocity discretization in (10)-(11) leads to a wrong limit as e goes to 0. Therefore, a change 
of variable in these velocity integrals is necessary before any discretization to ensure the right anomalous 
diffusion limit. 

In Case 1 (heavy-tailed equilibrium), we perform the change of variables w = e\k\v only in the second 
velocity bracket of (10) and do the same for (11). Then, a rectangle formula is applied to the obtained 
integrals. The first velocity bracket in (10) is computed directly using a rectangle formula, and the same 
is done for (11). We proceed similarly for Case 2 (degenerate collision frequency), except that we perform 
the change of variable w = e\k\v/v{y) at the same places. It is important to note that the discretizations 
in velocity are independent of e and we still have the first order (in time) uniform accuracy with respect 
to e; this statement is proved in [4]. 

In order to highlight the AP character of this scheme, we present in Fig. 2 the densities obtained with 
the initial data fo(x,v) = (l + sm(Trx))M(v),x £ [—1,1] at time T = 0.1 and with At = 10” 2 . In the case 
of the heavy-tailed equilibrium, we took /3 = 2.5 and in the case of the degenerate collision frequency, we 
took j3 = 0.5. In both cases, we consider u max = 5 and 200 points of discretization in velocity. 




Figure 2. The densities p v (T = 0.1, x) for the Duhamel schemes and the anomalous diffusion limit scheme. Left: the case of 
heavy-tailed equilibrium. Right: the case of degenerate collision frequency. 

Note that in both cases, the distribution function / can be easily recovered from the values of p u given 
by the previous scheme, using a Duhamel formulation of the kinetic equation for /. The same uniform 
accuracy property for / is inherited for that of p„. 
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